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Turbulent Torques on Protoplanets in a Dead Zone 

Jeffrey S. Oishi^'^, Mordecai-Mark Mac Low^, Kristen Menou^ 

ABSTRACT 

Migration of protoplanets in tfieir gaseous host disks may be largely responsi- 
ble for the observed orbital distribution of extrasolar planets. Recent simulations 
have shown that the magnetorotational turbulence thought to drive accretion in 
protoplanetary disks can affect migration by turning it into an orbital random 
walk. However, these simulations neglected the disk's ionization structure. Low 
ionization fraction near the midplane of the disk can decouple the magnetic field 
from the gas, forming a dead zone with reduced or no turbulence. Here, to under- 
stand the effect of dead zones on protoplanetary migration, we perform numerical 
simulations of a small region of a stratified disk with magnetorotational turbu- 
lence confined to thin active layers above and below the midplane. Turbulence in 
the active layers exerts decreased, but still measurable, gravitational torques on 
a protoplanet located at the disk midplane. We find a decrease of two orders of 
magnitude in the diffusion coefficient for dead zones with dead-to-active surface 
density ratios approaching realistic values in protoplanetary disks. This torque 
arises primarily from density fluctuations within a distance of one scale height 
of the protoplanet. Turbulent torques have correlation times of only ~ 0.3 or- 
bital periods and apparently time-stationary distributions. These properties are 
encouraging signs that stochastic methods can be used to determine the orbital 
evolution of populations of protoplanets under turbulent migration. Our results 
indicate that dead zones may be dynamically distinct regions for protoplanetary 
migration. 

Subject headings: accretion — MHD — planetary systems: formation — plane- 
tary systems: protoplanetary disks — turbulence 
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Introduction 



Detections of extrasolar planets now number in the hundreds ( iButler et aLll2006l ). Their 
orbital parameters reveal a nurnber of multi-pl anet systems in mean motion resonances 
( iLee fc Pealell2002l : |ji et al.ll2003l : iLee et al.ll2006l ) a nd many Jupiter a nd Neptune- mass ob- 
jects with major axes under an astronomical unit ( Gould et al.l 2006 ). These observations 
can be explained by orbital migrat ion arising from disk-planet interaction (ILin fc Papaloizou 
I979I : IColdreich fc Tremainelll98oh . 



Protoplanetary migration driven by gravitational interaction with a disk can occur in at 
least three ways. Type I migration occurs when protoplanets with mas ses < 10 — SOMm raise 
densi ty waves in the disk, which in turn exert torque on the protoplanet (iGoldreich fc Tremaine 
1980l ). This migration is usually inward due to the radial gradient s in pressure , tempera- 
ture, and sound speed from the global structure of typical disks (IWardI Il986l ). Type II 
migration is an essentially non-linear p rocess that occurs when a massive planet opens a 
gap in the surrounding accretion disk (IWardI 119971 ) . The planet then couples to the vis- 
cous evolution of the disk, migrating inward on the local disk accretion timescale. Type 
III migration occurs when material in the coorbital region of a roughly Saturnian mass ob- 
ject begins t o exert strong corotation t orque that scales with the migration rate, leading 
to runawav (IMasset fc Papaloizoul l2003l . though the existence of runaway is questioned by 
D'Aneelo et al.l J2005h on resolution grounds). We will here be solely concerned with Type 
I migration. 



[e.g.. 



For an ^ IMm object at 5 AU, th e time scale for Type I migration is r^jg ~ 8 x 10^ yr 
Tanaka. Takeuchi. fc Wmdll2002h . while protoplanetary disk s have lifetimes constrained 



to be less than tdisk ~ 10'' yr (e.g., Jayawardhana et al. 20od). This poses a significant 



challenge to the core-accretion scenario for gas giant planets ( iPoUack et al.l 119961 ). in which 
rocky cor es take ~ 10'' yr to gro w to ~ 10 — 30M® before rapidly accreting a gaseous 



envelope. iRice fc Armitagd (120031 ) note that if a core migrates stochastically, the timescale 
for core accretion drops by an order of m agnitude from this estimate. More recent studies 
( Hubickyj. Bodenheimer. fc Lissauil boosl ) show a core-accretion timescale of 1 — 5 x 10^ yr 
for Jupiter, still considerably longer than the migration timescale. 

Essential to this discussion of proto planetary migration is the disk itself. Disks around 
young stars are known to be accreting (IHartmann et al.lll998l ). Angular momentum must 
be transported outward to allow accretion. The strength of the transpo rt mechanism can be 
most simply parameterized by the dimensionless viscosity a = ut/ (csH) (IShakura fc Sunyaev 
I973I ). Models incorporating a viscosity characterized by a are known as a-disks. Nu- 
merical simulations of disk-planet interaction confirm the basic analytic predictions of the 
timescales for Type I migration in purely smooth, featureless, laminar disk models in- 
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eluding g-disks and the minimu m mass solar nebula (iMiyoshi et al.l Il999l : iMassetl 12002 
D'Angelo. Kiev, fc Henn"iii3 l2003h . 



However, protoplanetary accretion disks are far from featureless, laminar disks. Re- 
cently, two separate approaches to more complex modeling of planet-disk interaction have 

been identified and p ursued. First, more sophi sticated one-dimensional a-disk mo dels with 

radial opacity jumps jMenou fc Goodmanl[20oi ) and departures from isothermahty (jJang-Condell &: Sassek 



20051 ) show that Type I migration rates can be sigr iificantly altered by regions of lon g mi 



gration time that function as traps. Additionally, iPaardekooper fc Mellema I (l2006l ) have 
performed three-dimensional radiation hydrodynamics simulations of Type I migration and 
found that Type I migration can be halted and even reversed if the disk cannot radiate 
efficiently in the corotation region of the planet's orbit. Second, angular momentum trans- 
port most likely is actually driven by some f orm of turbulence, likely the saturated state of 



the magnetorotational instability (MRl) (see lBalbus fc Hawleylll998l . for a review). Under 



standing the effects of turbulence on migration requires explicit, three-dimensional, magne- 
tohydro dynamic (MHD) simulations. In a series of recent papers. Nelson, Papaloizou, and 



bedded within the turbulent flow ( 


Papaloizou & Nelson 


200c 
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2003; 


Papaloizou. Nelson. &; Snellgrove 


2004 
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These models revealed that for an object with 



< 3OM0, torque fluctuations from turbulent density perturbations cause the protoplanet to 
follow a random walk, possibly significantly lengthening its lifetime. 



Johnson. Goodman, fc Menoul (120061 . hereafter JGM06) have argued, using semi-analytic 
models, that these fiuctuations will decrease the lifetimes of most planets, but allow a small 
number of them to scatter to large radii and thus survive. Their approach is based on a 
Fokker-Planck formalism that treats Type I migration as advection and the orbital motions 
induced by stochastic torques from turbulence as an independent diffusion process. Such a 
model depends on details of the turbulence that can only be provided by numerical simula- 
tion. 

Because protoplanetary disks are cold and dense at their midplanes, the ionization 
fraction is a strong function of height. The MRI is stabilized by Ohmic diffusion for magnetic 
Reynolds numbers Reu ^ 10^ for zero-net flux poloidal fields, though the actual value 
depends strongly on the geometry of the fleld (IFleming. Stone, fc Hawlevll2000l) and possibly 
the presence of artiflcial resistivity (IBrandenburg et al.ll2004j ). T herefore, the M RI can not 
operate at all heights at radii where planet formation likely occurs (Gammie 1996h. Modeling 



of the ionization structure of the disk leads to a three layer description (I Gammid Il996 



Sano et al.ll2000l : lFromang et al.ll2002l : ISemenov et al.ll2004J : llkner fc Nelsonll2006f ) with MRI 
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turbulence occurring in active zones above and below a dead zone at the midplane. 

Our purpose in this paper is twofold. First, we wish to understand the effect of a 
dead zone on the migration of a low-mass protoplanet at the midplane. Second, we wish to 
further quantify the torques on a protoplanet caused by MRI turbulence, with and without 
dead zones, in order to verify the assumptions of and provide better parameters for the 
Fokker-Planck model of JGM06. 

§ [2] describes our models, § [3] gives a brief overview of the turbulence and dead zones in 
the plasma, §|1] provides analysis of the torques on protoplanets in real and Fourier space, §[5] 
discusses the applicability of our results to real disks and statistical treatments of ensembles 
of migrating protoplanets, and § [6] presents our conclusions. 



Models 



Torques produced by turbulent density pert urbation s comp ete with Type I migration 
in the turbulent migration scenario proposed by iNelsonI (120051 ). In order to focus on the 
characteristics of the turbulence, we exclude active Type I migration in our models by 
considering a zero-mass test particle fixed at the center of the box. A test particle does 
not raise density waves in the disk, allowing us to isolate the torques on the planet from disk 
turbulence. We work in a local frame, allowing for maximum resolution of the disk near the 
planet — where, as we will show, most of the turbulent torque originates. Our test particle 
approach has the additional advantage of allowing us to study turbulent migration in the local 
frame without worrying ab out spurious density wakes from the (shearing) peri odic boundary 
conditions we employ (eg. iPapaloizou et al.ll2004j : iNelson fc Papaloizoul 120041 ). However, in 
doing so, we implicitly assume that Type I and tu rbulent torque s are separable and additive. 
It is not clear that this assumption is valid, as iNelsonI (j2005( ) found that explicit Type I 
migration in turbulent models differed from the linear combination of stochastic migration on 
massless protoplanets and a Type I migration rate from laminar disk calculations, suggesting 
that Type I migration is itself modified by MRI turbulence. 



2.1. Numerical Method 



We use the Pencil Codj^, a spatially sixth- order and temporally third-order finite differ- 
ence MHD code (IBrandenburg &: Doblerll2002l ) designed to study weakly supersonic turbu- 



available at |http://www. nordita.dk/software/pencil-code/ 
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lent flows. The Pencil Code solves the MHD equations in non-conservative form, monitoring 
conservation as a check on the quality of the solutions. 



We work in the shearing box approximation (e.g.. lGoldreich &: Lynden-Belllll965l : lHawley. Gammie. fc 



19951 ). in which one considers a small Cartesian box around a fiducial point in the disk at 
radius R, and expands the gravitational potential from the central protostar to leading order 
in H/r. The coordinate frame is oriented such that x is the local radial direction, y points in 
the azimuthal direction, and z is mutually orthogonal. We solve for the departures u from 
the mean Keplerian shear flow, u^^ = —3/2Qx. The shearing box MHD equations are thus 

5,p + «f9,p + V-(pu) = /, (1) 

VP JxB „ _„ . nu. 



dtu + u ■ Vu + u^'^^dyU 



P 



P 



fu + 2Quyyi - 



-y-n'2z + C(VV-u) (2) 



dtA + uf^dyA 



SUA 



y , 



-x + u X B + r7(2;)poJ + fA (3) 

where f2 is the rotation rate of the box, /p, fu, are stabilizing hyperdiffusivities (see 
below), A is the magnetic vector potential, J = V x B/po is the current density, rj{z) is a 
height-dependent resistivity, and the other symbols have their standard meaning. The extra 
term u^y^ dy in all dynamical equations is a result of subtracting the Keplerian shear from the 
velocity, as is the {?)/2)VLAy term in equation ([3]). Finally, we close the set with an isothermal 
equation of state, 

P = cIp, (4) 

where is the isothermal sound speed. 



Because the Pencil Code uses spatially centered finite differences, there is no formal 
diffusion in the algorithm, so explicit diffusion operators must be included in all dynamical 
equations in order to ensure numerical stability. In this work, we use sixth order hyper- 
diffusion for continuity, momentum, and induction, in addition to physical Laplacian (i.e., 
second order) diffusion. Sixth order hyperdiffusion has the form i/gV^u and damps signals 
only near the Nyquist wavenumber at the smallest scales of the box, allowing larger modes 
to remain largely unaffected. Our form ulation conserves raass an d momentum by construc- 
tion. For details, we refer the reader to lJohansen fc Klahrl (120051 ). We capture shocks using 
a von Neumann artificial viscosity implemented as a bulk diffusion with variable coefficient 
= c(max[(— V • u)-i-]) (where the max is taken over 3 zones) on the continuity and momen- 
tum equations (IHaugen. Brandenburg. &: Meell2004l ). 



For our dead zone runs, we hav e modified the Pencil Co de to include a height-dependent 
resistivity with the profile given by [Fleming fc Stond (120031 ): 

f z^\ f I 
?7o exp — - exp 



ri{z) 



^CR 



2V^ 



dz 



(5) 
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where T^cr — 100 g cm~^ is the stopping depth of cosmic rays and Sq is the surface density 
of the disk. We adopt a value for Sq/Scr = 30, consistent with FS03 and roughly twice 
the minimum mass solar nebula value at 1 AU. Note that FS03 incorrectly identify rjo as 
the midplane resistivity, rather than rj^id = ''7o exp (Eo/4Ec:r), though this error does not 
carry through their calculations of the magnetic Reynolds number (J. Stone 2006, private 
communication) . 



2.2. Initial and Boundary Conditions 



We initialize the disk in hydrostatic equilibrium in the vertical direction, leading to 
a density profile p{z) = poexp{—z'^/H'^), where H = -\/2cf7^ is the scale height of the 
disk. The initial magnetic field strength is parameterized by the maximum plasma (3 = 
2iiQPg/ Bq = 400. The field geometry is a zero-net flux vertical field Bz{x) = Bq sm{2nx/Lx). 
We measure length in units of the disk scale height, H, time in units of (although we 
plot quantities in terms of the orbital period torb = ^n/Q), and set the vacuum permeability 
Po = 1. Our computational domain has size (a; x ?/ x 2;) = 1 x 4 x 4, with resolutions for the 
models given in Table [H 

In these units, all simulations have a sound speed Cg = 7.071 x 10~^, rotation rate 
Q = 10^^, midplane density po = 1? and random perturbations in u with amplitude 1 x 10^® 
to seed the MRI. In order to vary the size of the dead zone, we vary 770 in equation [5l which 
in turn varies the midplane resistivity. We parameterize the size of the dead zone with the 
magnetic Reynolds number at the midplane Rcm = / (Vmid^) ■ 



Our boundary conditions are periodic in y, shear periodic in x (iHawley et al.lll995l ). 
and periodic in z. While the last is not rigorously accurate for a stratified disk such as 
ours, previous work (IStone et al.lll996f) as well a s our own direct comparison to vertical field 
boundary conditions (IBrandenburg et al.lll995l ) suggest that it does not make a significant 
difference for the density perturbations driven by the MRI. 



2.3. Validity of Zero-mass test particles 

We monitor the gravitational force on a zero-mass test particle fixed at the center of 
the grid. This can be directly translated into a torque for a particle that does not move 
significantly in radius over the course of the simulation. While fixing the position of the 



planet over the ~ lOOtorb duration of the simulation is an approximation. Figure 6 of iNelson 



(120051 ) shows that zero-mass objects vary their semi-major axes by no more than about 10% 
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over a similar time span. 



If we make the assumption that Type I and turbulent migration are separable, then our 
analysis will apply more generally to objects of finite mass. However, even if these effects 
do couple, there exist a class of objects for which our results are valid. Here, we derive a 
range of masses for objects large enough to neglect gas drag and still small enough to cause 
negligible density waves. 

Stokes' drag law states that the stopping time of an object larger than the mean free 
path of the gas it is moving through is 



2psa 



9/i 



(6) 



where ps and a , are the density and r adius of the solid object and /i is the dynamic viscosity 
of the gas (e.g. IWeidenschillingl 119771 ) . 



In order to determine the mass range in which an object can have both drag force 
and gravitational back reaction on the disk neglected (ie, a massless particle on a fixed 
Keplerian orbit), we consider the range of sizes a body with solid density roughly that of 
Earth, ~ 3 g cm~^, will have both drag times and Type I migration times on the order 
of the lifetime of a protoplanetary disk, taken here to be 1 0^ yr. The dynamic viscosity 
of a pure hydrogen gas is /i = 5.7 x 10~ ^r~^/^g cm~^ s~ ^ ( Land Il980l ). Using an a-disk 



10^2 and M = 10"^ M© yr" 



Ilgner fc Nelson! (120061 ) give a midplane value of 
100 K. The lower limit, Tf > 10^ yr gives a radius ~ 5.1 x 10'^ cm which for a 



with a 
T 

spherical object of the above density yields m > 4 x 10^^ g. Type I migration times scale 
inversely with mass, rj ~ 8 x lO^(m/M0)~^ yr (ITanaka et al.ll2002l ). and so the upper mass 
limit is rrimax ^ 10^^ g. This is a broad range of protoplanet masses, 10^'* g ^ ""^ ^ 10^^ g, 
and suggests our results are relevant even if Type I and turbulent migration are coupled for 
more massive objects. 



3. Plasma Dynamics and Dead Zone Parameters 



In our models with finite Rcm, magnetic energy rapidly grows near the midplane after a 
few orbits , due to the coope r ation of magnetic reconnection and shear, in agreement with the 
results of [Fleming fc Stond (120031 ). By 10 orbits, MRI turbulence sets in above and below 
the midplane, forming active layers surrounding a quiescent (but not motionless) dead zone 
that persists to the end of the simulation (Figure [1]). We define the dead zone to be the 
region where the horizontally averaged Maxwell stress —{BxBy)xy, has fallen by an order of 
magnitude from its average value at ±1.5//. (We use unadorned angle brackets to denote 
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volume averages, while spatial averages over lower dimension and temporal averages have 
identifying subscripts.) Using 1.5H as our baseline, rather than the grid boundary at 2H, 
avoids contamination by the slightly elevated value of the Maxwell stress that we see at the 
boundary. We have confirmed this effect is a result of the periodic boundary conditions 
adopted with simulations using other boundary conditions. We calculate the total surface 
density in the active zones and compare that with the surface density of the dead zone 
T,£) in Table [H 

In a shearing box, symmetry restricts the net transfer of angular momentum, and thus 
the accretion rate, to zero. In order to measure the ability of the turbulence to drive accretion, 
we measure the Maxwell stress, as well as the Reynolds stress {pu^Uy). Th ese quantities 



measure the radial transport of angular momentum by the turbulence (e.g. iBrandenburg 



19981 ). Figure [2] shows the turbulent stresses normalized by the initial midplane pressure 
for each of the medium resolution runs. We see Reynolds stresses significantly exceeding 
Maxwell stresses in the midplane for all of our dead zone models. FS03 demonstrated that 
these correlated motions are caused by momentum flux from turbulent overshoot at the 
boundary layer between dead and active regions. 

Our models confirm the finding of FS03 that a significant Reynolds stress remains within 
the dead zone, even as the size and depth of the dead zone, as defined by the dropping 
Maxwell stress, increase. Because the criteria used by FS03 to define the dead zone width 
are not entirely clear, it is difficult to make detailed comparisons between our results and 
theirs. However, applying our definition of a dead zone to Figure 3 in FS03, which shows 
Maxwell and Reynolds stresses for a model with Rcm = 100, the dead zone extends from 
\z\ < 0.5. In our comparable model, the dead zone is located between —0.35 < z < 0.38. 
Given that FS03 report that their Rcm = 10 model was unable to sustain MRI turbulence 
while we were able to run a self-sustaining Rcm = 3 model, it is possible that the difference 
in dead zone widths is a result of the different dissipative properties of our respective codes. 

The Rcm = 3 case demands special attention. Its vertical distribution of stresses is 
notably different than those of the other non-ideal runs. First, the Reynolds stress at the 
midplane is substantially less than in the Rcm = 30 and 100 cases, which differ by a factor 
of a few. Second, the Maxwell stress appears to reach a minimum below 10^^ and remain 
constant over 1^1 < 0.5, though it occasionally becomes slightly negative and is excluded 
on the log plot. This behavior occurs in both 32R3 and 64R3 mod els, as well as for a run 



in double precision. We note that Figure 5 of [Turner et al.l (120061 ) shows similar behavior 



(though not the negative —B^By) in dead zone models computed with the Zeus code. 
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3.1. Resolution Study 

In order to study the effects of numerical resolution on our results, we ran our Rcm = 30 
model at three resolutions, and all others at two. An overall picture of the MRI is given by 
its saturation energy or average viscous a value, both of which are shown in Figure [31 A 
more relevant quantity for our purposes is the temporal root-mean-square (RMS) azimuthal 
gravitational force ((-F^)*)"'^ (i.e., torque; see §|4]). This quantity is directly responsible for 
orbital migration and is given in Table [H Figure [3^ shows the RMS torque and dead zone 
column density ratio Sa/S^ as a function of resolution. While both are resolution dependent, 
it appears that their values are converging between the highest two resolution runs. The 
dimensionless effective viscosity a = {{pu^Uy) — {Br^By)) / Pq and maximum Mach number 
are given in the upper two panels of Figure [3)d. Both quantities refle ct the well-known result 



that higher resolution MRI runs lead to lower turbulent stresses (IFromang fc Papaloizou 



20071 ) 



The magnetic energy for the high resolutio n run drops by nea rly a factor of 7 between 



orbits 30 and 50, similar to the simulations of IStone et al.l (119961 ). Decreases in magnetic 
energy at similar times and on similar timescales occur at all resolutions, although not 
with such large magnitudes. In the lower resolution models run to longer times, these 
drops eventually stop occurring, typically by orbit 100 and the model maintains a relatively 
constant magnetic energy level. The timescale for such an adjustment is curiously long; 
we hope to examine it more closely in future work. For our purposes here, it is sufficient 
to note that while resolution effects are present, the du ration of the simulation may be a 



greater factor in obtaining robust statistical results (see [Winters. Balbus. &: Hawleyl 12003 
for a discussion). 



4. Migration Torques 

Turbulent overdensities in a protoplanetary disk exert torques on a protoplanet at radius 
R with strength 

r = rfJ/dt = Rx F. (7) 
The sum of gravitational forces from all turbulent density perturbations is 

where ri is the distance to each gas zone with density pj, = AxAyAz is the zone 
volume and the sum is taken over all zones in the computational domain. In particular. 
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orbital migration is caused by a change in angular momentum J, which for a circul a r orbi t 



in a spherically symmetric potential is simply J^. Following iNelson &: Papaloizoul (120041 ). 
we consider only a protoplanet at constant R — the center of our shearing box, and thus 
dJz/dt = oc Fy. In the following analysis, we consider only the properties of the scalar 
Fy and thus only track migration on circular orbits. We scale Fy in units of 27rGS, which 
is the force per unit mass felt by a particle suspended a small distance above the center 
of a disk with constant E. We sample this quantity every 100 timesteps throughout the 
calculation. While the timestep is dynamically determined and thus fluctuates over the 
simulation, it maintains an average value of ~ 10~'^torb with an RMS at least an order of 
magnitude smaller, giving a sampling rate of roughly At ~ lO^^torb- 

Figure H] shows Fy as a function of time for simulations with varying values of -Re a/. 
The RMS value of the force declines with increasing dead zone size. This decrease occurs 
because it is the MRI in the active layers that drives the density fluctuations, which are in 
turn responsible for the azimuthal force. For the remainder of the analysis, we consider only 
torques for times t > 25toih in order to avoid any of the transient features from the onset 
of the MRI evident in Figure HI We will refer to the value of Fy as "torque," although it is 
formally a force per unit mass. 

Two critical assumptions of the Fokker-Planck model for planetary migration of JGM06 
are that the turbulent torques have temporal stationarity and finite correlation time. We 
test each of these assumptions for the torque distribution in our simulations. 



4.1. Torque Distribution 

In order to determine time stationarity, we need the distribution of torques over various 
fixed time intervals. We can then ask if the samples in each interval are consistent with 
being drawn from the same distribution as the other intervals. We separate the torque time 
series into seven lOtorb blocks from 25 — lOOtorb- Figure shows that the mean of each 
interval remains near zero over time, and the RMS remains nearly constant. JGM06 require 
that the torque distribution be such that ST(t, jy has no time dependence, and this plot 
demonstrates this for our study. However, the means over each lOtorb bin vary about zero 
with amplitudes a factor of < 0.2 times the standard deviation over the entire 75torb sample, 
suggesting that the underlying distribution of torques does vary slightly over this short an 
interval. This does not affect the stochastic migration presented here, as such an interval is 
much shorter than a typical diffusion time in the JGM model. 

Additionally, we are interested in the particular distribution of these torques: the central 
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limit theorem assures us that the cumulative effect of any random process limits to a Gaussian 
distribution, but can we expect to have a coherent Gaussian over a given (short) period, given 
the unknown distribution of turbulent forces? Figure E] shows the distributions of torque over 
10 orbit periods for each of the dead zone widths. Each period represents ~ lO'^ samples. 
There is a significant change in the distribution for the Rcm = 3 case, although it is consistent 
with the trend of decreasing standard deviation with decreasing Rcm- Qualitatively, even 
over such short periods, the torques appear Gaussian distributed, in the sense that they 
are centrally concentrated on zero and roughly symmetric (figure E]) . We quantify this by 
computing the skew 



S 



1 ^ 
N ^ 



a 



-I 3 



and kurtosis 



K 



1 ^ 
N ^ 



(9) 



(10) 



where a is the standard deviation of the samples, given in tables [2] and [3] respectively. A 
normal distribution has skew and kurtosis both equal to zero, while an exponential distribu- 
tion has kurtosis of three. These tables show that over 10 orbits, the distribution varies from 
normal signi ficantly — the varia nce from normal is given using the standard a^ar = 
for kurtosis (jPress et al.lll992l ). This is directly related to the chaotic nature of the MRI: 
small deviations in initial conditions pr esent at the start of each lOtorb block will cause sig- 
nificant deviations (IWinters et al.ll2003l ). which we observe as differences in distribution. We 
note here that this analysis does not affect the JGM assumptions-their model does not re- 
quire a Gaussian distribution (shown above)-it merely underscores the short-term properties 
of the MRI. 



4.2. Fourier Analysis 



Figure [7] shows the temporal power spectra of the torque for runs with varying dead 
zone widths. Increasing dead zone width has little effect on the shape of the spectrum, 
changing only the total power. This suggests that the turbulence produces a similar spec- 
trum of density fluctuations, regardless of the width of the dead zone, as we wou ld expect 



since the underlying instability is not significantly altered by Ohmic resistivity (jjini Il996 
Fleming et aPboOOT l. 



The low frequency power reported by iNelsonI (120051 ) is absent in our models, done in 
the local frame. This suggests that whatever imparts a long-term memory to the turbulence 
in his simulations is related to the global structure. This conclusion was also noted in 
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Nelson fc Papaloizoij (120041 ) on the basis of running time averages for the migration torque 
in their local simulations. They find that running averages of the torques tend to zero, while 
for low mass planets in global simulations, they do not. Our stratified but non- resistive 
models show a similar convergence, which tends to smaller values at late times for increasing 
dead zone width (Fig. [8]). 

As noted in § HI our time series data is not evenly spaced in time. To correct for this, 
we interpolate to an even time grid for Fouri er analysis. W hile this can be dangerous for 
astronomical time series with large gaps (eg.. IScargld 119821 ). we have confirmed that in our 



case, interpolation is essentially identical to a Lomb Normalized Periodogram (jPress et al. 



19921 ). which produces power spectra (but not correlations) for unevenly sampled time series 
data. 

The existence of discrete peaks in the power spectrum indicates the presence of char- 
acteristic timescales for the torquing of the protoplanet by MRI turbulence. A more direct 
measure of this timescale can be found with the autocorrelation function, which we calculate 
by taking the inverse Fourier transform of the power spectrum. Figure M shows autocorre- 
lation functions for each of the medium resolution runs. There is almost no correlation for 
lags longer than about a dynamical time. Thus, in a local frame at least, turbulent torques 
are correlated only for a short time, typically ~ O.Storb (see Tabled]). 



4.3. Locality 

It is interesting to understand where the forces relevant to turbulent migration come 
from. Are they dominated by local forces, or do increasingly large perturbations from large 
scale motions dominate the total force? Figure [10] shows the total instantaneous force con- 
tributed by concentric spherical shells centered on the protoplanet at the box center, cal- 
culated as AFy{r) = J^i^^yixijiii, Zi), for < xj + yf + zf < (r + dr)^, at four random 
times for each model. All panels show a sharp drop-off outside of r ~ H, validating the 
local approach and demonstrating that nearby perturbations tend to dominate over larger 
turbulent structures. 

However, the shearing box imposes a cutoff on the outer scale of MRI turbulence that 
precludes very large eddies from forming. In particular, at scales large enough for the disk to 
act as a two-dimensional flow, eddy formation may occur from the inverse cascade of energy 
expected in two-dimensional flows. We cannot rule out the possibility that such structures 
may affect protoplanetary migration. JGM06 argue that large eddies will dominate the 
diffusion coefficient in a turbulent flow because they have longer correlation times. Our 
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simulations do not address this, as we do not have information about the correlation time as 
a function of scale. However, we have run a model with twice the radial (x) extent that shows 
a similar drop-off at a scale of H ensuring that our result is not simply a numerical effect 
of the small radial size of our box. In addition, we have some insight from the observation 
that the correlation times differ very little through the sequence of increasing dead zone sizes 
(see § 14. 2p . This suggests that the total effect of the various size density structures does not 
change, even as the various eddy sizes change as a result of the dead zone. Future work will 
address this point. 

Overdensities at the midplane dominate the torque not only in the fully ionized model, 
but even in models with large dead zones, for two reasons: obviously, they are closer to the 
protoplanet, and secondly, less obviously, the higher total density at the midplane means that 
even small perturbations there produce absolute overdensities higher than in the turbulent 
active layers. Figure [TT] shows images of the turbulent overdensities, p — po{z) at t = lOOtorb 
for all values of Rcm- Subtracting off the initial exponential profile po{z) simply removes 
the large symmetric component due to the unperturbed density distribution, allowing the 
density fluctuations that contribute to the torque to be more easily visualized. The scale 
bars show the general decline in density perturbation amplitude. Two features are worth 
noting: the increasing dominance of wave-like structures in these torque images, and the fact 
that the torque drop-off does not come from the MRI moving further away from the planet. 
While that does happen, the strongest density perturbations remain near the midplane, 
regardless of the presence of the MRI at that location. The decrease in torque strength for 
lower magnetic Reynolds numbers is due to the overall lessening of turbulent energy because 
of both thinner active layers, and reduced energy at saturation due to the higher Ohmic 
resistivity present even in the active layers. For our ideal MHD runs, 6p is consistent with 
the relation 

{Spr\.{sPLS" 



( ISano et al.ll2004l ). suggesting that their result also applies to stratified disks. For our isother- 



mal runs,(P) = (c^p) = (p) and the relation reduces to (5p^)°'^ ~ (SP^ag)^ ^ I'^s- Table H] 
shows the relationship between volume averaged density and magnetic pressure perturba- 
tions for our medium resolution models. As the dead zone size increases so too do the density 
perturbations relative to the variance of magnetic pressure. This is consistent with the ob- 
servation of density waves in the dead zone: the active zones drive density perturbations in 
regions with negligible magnetic pressure, and so volume averages should show just such a 
trend to higher ratios of the two quantities. 

The wave-like structures that we see in th e simulations may be the lo^v-m sp iral density 



waves reported in earlier shearing box studies (iHawley. Gammie. fc Balbuslll996l ). The wave 
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pattern appears similar in all cases because it is driven by the MRI. We switched off the 
Lorentz force in one of our Rcm = 30 simulations, and allowed the evolution to continue for 
about 6 orbits. With no MRI driving in the active layers, the turbulence quickly dissipates. 
However, during the decay, the wave pattern remains, although its amplitude drops by 
roughly a factor of 8 in those 6 orbits. This suggests that the wave pattern may be a natural 
resonant state of either real disks or merely the shearing box. Although the images reveal 
wave-like structure, the temporal power spectra of torque do not differ significantly between 
fully turbulent and dead zone models. Simulation s of stochastic excitation of solar p-modes 
(INordlund fc Steinll200ll : IStein fc Nordlundl l200ll ) tend to show strong temporal signatures 
of the excited waves, while our power spectra do not show any obvious wave signatures. We 
will investigate these structures further in future work. 



4.4. Parameterization for Stochastic Models 

One of our objectives in this work is to provide parameters necessary for refining the 
statistical treatment of migration by JGM06 and to understand the predictions it makes 
when we include dead zones. The extreme computational cost of three-dimensional MHD 
simulations ensures that the fully global simulation including several scale heights and a few 
tens of AU in radius evolved to the ~ lO^torb required to follow planetary migration directly is 
unlikely in the foreseeable future. Thus, a promising tool in the study of planetary migration 
is a stochastic method in which one considers the fates of a population of planets subject to 
random torques from the disk turbulence. JGM06 have constructed an advection-diffusion 
equation from the Fokker-Planck equation in which standard Type I migration is modeled 
as advection and turbulent torques are modeled as diffusion (note the important correction 
in the erratum). Understanding the long term evolution of a population of planets in this 
model requires a measurement of the diffusion coefficient due to turbulent torques in regions 
of disks with varying resistivity profiles. 

JGM06 derive this diffusion coefficient in terms of the independent variable J, the orbital 
angular momentum of a protoplanet at a given radius Vp, rather than the radius itself. From 
a standard Fokker-Planck derivation, JGM06 define the diffusion coefficient as 

/oo 
6T{t - r/2, J)6T{t + r/2, J)dT ~ 6T{t, Jfr^, (12) 
-oo 

where Tc is the correlation time, estimated here from the autocorrelation function (§ 14. 2p . and 
the overbars represent ensemble averages. Direct calculation of the integral from our time 
series for Fy oc r(t, J) does not yield satisfactory results, probably because each simulation is 
only one realization of the torque fluctuations, and D{J) is well-defined only for an ensemble 
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of such realizations. In lieu of direct calculation, we make use of the rough equivalence 
D{J) ~ {ST{t, JY)tTci where we have substituted a time average for the ensemble average in 
the original. 

JGM06 collect all the uncertain properties of turbulence that need to be derived from 
a numerical model into a single parameter 



torb V 0-046 ) (0.046)%Mj 
where the dimensionless coefficient 



Co ^ (14) 



We calculate Cd using our simulated values of Fy, assuming fixed r^, and allowing and Mp 
to cancel between 6r oc VpMpFy and the denominator. Figure [12] shows e as a function of dead 
zone size in our various simulations. We assume that the ratio of surface densities S/j/S^ 
mainly determines the strength of fiuctuations in Fy. T his facilitates a comparison with 



the s i zes of dead zones pred icted by ionization models (e.g. ISano et al.l 120001 : iFromang et al. 



2OO2I : lllgner fc Nelsonll2006h . 



Figure [12] includes a resolution study. The circled values show the same model at all three 
different resolutions, while the two sets of models with extreme values of -Re a/ were run at the 
lower two resolutions. (Remember that the same value of Rcm produces different dead zone 
thicknesses at different resolutions, as shown in Figure [3](^aj.) All of these comparisons show 
that our derived diffusion parameter e drops steadily with increasing resolution, suggesting 
it is less well converged than the average torque or dead zone thickness (Fig. [^('aj). Thus, 
our estimates of the diffusion coefficient and related quantities appear likely to be firm upper 
limits to the true values. 



JGM06 derive a fiducial value of e ~ 0.5 from the ideal MHD simulations of Nelson 



(I2OO5I ). Our local ideal MHD simulation, on the other hand, gives a value of e ~ 10~^ at 
medium resolution. On the other hand, JGM06 derive a fiducial correlation time of ~ 0.5 
from Nelson's simulations, which is quite close to our ~ 0.3 value. This suggests that it is 
differences in the torque magnitudes, rather than correlation times, that drives the difference. 
The torque magnitudes are driven directly by the MRI perturbations, which in turn may 
differ from N05 due to his lack of stratification, by the difference in MRI strengths (his 
models have a ~ 5 x 10"^, larger than ours by a factor of a few), or by a difference in 
density perturbations between global and local models. In our models with dead zones, we 
find 10"^ < e < 10-2 over 3 < i?eAf < 00. 
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In the disk model that we employ, even our ideal MHD model gives a value of e that 
suggests a terrestrial-mass planet in even our most turbulent model will be in the advection 
dominated regime, as shown in Figure 7 of JGM06. (Note that the fiat parts of the curve 
in that Figure correspond to advection.) However, even this reduced value of e is sufficient 
to drive diffusive migration for protoplanets in the sub-Earth mass range. For such planets, 
dead zones may shift migration from diffusion to advection as they pass radially in from a 
fully active region with stronger turbulence. 



4.5. Dead Zone Sizes 



Given the decrease of e( J) with increasing dead zone size that we find at all resolutions, it 
is reasonable to ask how close our models come to a realistic dead zone thickness. Previous 
studies of the ionization fraction in disks robustly show that the presence of small dust 
grains (< 1 /im) with interstellar abundance will cau se dead zones to be orders of magnitude 



larger in column den sity than their active zones (iSano et al.l l2000l : iSemenov et al.l 12004 
Ilgner fc Nelson|[2006l l Our i?eM = 3 models only reach a surface density ratio S/j/S^ — 10, 



making dynamical models with significant submicron dust grains well outside of current 
computational limits. However, once these dust grains either sediment to the disk midplane 
or coagulate, ionization occurs much more rapidly and the dead zone size shrinks. This may 
occur in a small fraction of the disk lifetime, although a population o f new small grains with 
differ ent properties could then be produced by collisional processes ( jPullemond fc Dominik 
20051 ). For an a-disk with M = 10~^ Mp^/ yr and a = 10~^, the surface density ratios become 



Eo/Ea = 9.1 at 1 AU and 1.1 at 5 AU fllkner fc NelsonI 120061 ) . While this model assumes 
an a-value somewhat greater than we find, it suggests that simulations such as ours may not 
be so far from realistic dead zone sizes. Our results indicate that this ratio infiuences the 
magnitude of turbulent torques and thus could determine the transition between advection- 
and diffusion-dominated migration regimes for a variety of protoplanetary masses. 



5. Discussion 

Our results indicate that in the thick dead zone that is expected around 1-5 AU, tur- 
bulent torques will be reduced by at least two orders of magnitude from those expected in 
a fully turbulent disk at the same radii, allowing a uniform inward Type I migration mode 
to become dominant even for protoplanets with Mp <C M^. A potential consequence of this 
is that if such protoplanets form further out in the disk where the dead zone is thinner, 
the inner, thick dead zone region of the disk could act as a sink for protoplanets following 
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random walks driven by the more vigorous turbulence expected further out. 

Because our simulations are local, we can only provide ST{t, J) at one value of angular 
momentum J (that is, one radial location in the disk). A radial gradient of the diffusion 
coefficient D{J) contributes to advection and does not affect our main results on dead zone 
diffusive migration. Indeed, a strong gradient in the diffusion coefficient might be present 
at the radial boundary between an active and dead zone in the midplane. Most models of 
protoplanetary disk ionization suggest that there will be two such boundaries, one in the 
outer disk where the column density increases to the point that protostellar X-rays can no 
longer penetrate to the midplane, and one in the inner disk where thermal ionization begins 
to become important. JGM06 note that the orbital flux Fj of the planet distribution function 
/ can be written as 

Fj = (r-djD)f-Ddjf, (15) 

emphasizing the contribution of djD to the effective mean migration. In a region where 
protoplanets cross from an active region to a dead zone in the midplane, e, and thus D 
(eq. [T3l) changes by several orders of magnitude (Fig. [T2|) over a region of uncertain but 
likely sub-AU width. At the inner boundary of the dead zone, djD < contributing to 
outward mean migration, while at the outer boundary, djD > 0, contributing to inward 
mean migration. Thus, the random walk migration induced by a fully turbulent disk tends 
to make the dead zone act as a trap for sub-terrestrial mass protoplanets. 

The magnitude of this effect is difficult to estimate from our simulations because we 
ignore background gradients in the disk, and thus our measurements of dead zone and fully 
active turbulent torques are at equivalent radial positions in the disk. Nonetheless, future 
global and statistical models should be able to predict if the effect overcomes Type I migration 
in these regions strongly enough for the dead zone to function as an efficient trap. 



Recently, [Turner. Sano. &: DziourkevitchI (120061 ) have performed simulations including 



a dynamic r]{x, y, z, t) in which the dead zone disappears, at least from the point of view of 
angular momentum transport: their results suggest that while the MRI remains suppressed 
below \z\ ~ 0.5H, there is still significant Maxwell stress from mean fields that build up due 
to the vertical expulsion of flux and the regeneration of azimut hal field f r om sh ear. They 



suggest that these fields may grow large enough to invoke the iTerquemI (120031 ) migration 
mechanism in which magnetic resonances close to the planet may reverse Type I migration. 
However, these results — while provocative — should not significantly alter our main results. 
In our case, we isolate the turbulent contribution to migration, ignoring the gravitational 
back-r eaction of a pro toplanet on the disk. Any effects of the midplane mean field, such 



as the iTerquenj (120031 ) mechanism, are the refore absent b y cori struction. Indeed, since the 



MRI is suppressed at the midplane in the [Turner et al.l (120061 ) simulations, the turbulent 
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overdensities will be reduced, as in our simulations. 

Our work depends on two major approximations: the restriction of our simulation to 
the local, shearing box, and the adaptation of zero-mass test particles fixed at the center 
of the box. Because of the latter assumption, the possibility exists that our results un- 
derestimate correlations because a moving particle would accelerate in the direction of the 
largest current density fluctuation, thus increasing the torque. However, the density struc- 
tures will still have lifetimes on the order of an orbit and they will not have different spatial 
distributions. This means that in the lifetime of the fluctuation, the protoplanet will not 
have time to move significantly from its current position. Estimating the distance by an 
impulse approximation, the distance moved over the lifetime of a perturbation is roughly 
x/H 10~^-a trivial amount. Any global dynamical structures associated with the MRI 
are absent from our models, including large scale eddies that could change our results if they 
are sufficiently massive. Large scale dynamical features may also affect the correlation times 
of torques. Despite these caveats, we generally expect our results to be qualitatively valid in 
demonstrating that dead zones are regions of significantly reduced diffusive migration. Dead 
zones may also have other co nsequ ences for planetary migr ation such as those considered by 



Matsumura fc Pudrita (120031 ) and iMatsumura et al.l (120071 ) . 



6. Conclusions 

We find that a dead zone in a protoplanetary disk can reduce the magnitude of turbulent 
torque fluctuations on migrating protoplanets by at least two orders of magnitude from a 
comparable fully ionized region. The thick dead zone expected at 1 AU, assuming a dust- 
depleted disk, can have significant consequences for migration of planets with Mp <^ 
driven by turbulent density fiuctuations. 

The central assumptions of the stochastic diffusion model of protoplanetary migration 
of JGM06 amenable to test in local simulations seem well-founded. The torques have a 
finite correlation time and have stationary distributions. These fiuctuations maintain nearly 
constant correlation times, independent of their amplitude. The stochastic diffusion model 
suggests that if protoplanets form outside of the dead zone and migrate in a diffusive manner, 
the dead zone may act as a sink. 
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Table 1. Basic parameters for all runs. 



Name 


Resolution 




t max / torh 








32Rinf 


32 X 128^ 


oo 


100 




5.17(-3) 


0.31 


32R100 


32 X 128^ 


100 


100 


1.65 


2.63(-3) 


0.31 


32R30 


32 X 1282 


30 


200 


1.13 


2.27(-3) 


0.39 


32R3 


32 X 128^ 


3 


100 


0.177 


5.17(-4) 


0.24 


64Rinf 


64 X 256^ 


oo 


100 




3.23(-3) 


0.31 


64R100 


64 X 256^ 


100 


100 


1.37 


1.90(-3) 


0.43 


64R30 


64 X 2562 


30 


100 


0.860 


1.49(-3) 


0.29 


64R3 


64 X 256^ 


3 


100 


0.102 


3.45(-4) 


0.33 


128R30 


128 X 5122 


30 


50 


0.660 


1.29(-3) 


0.25 



'^Magnetic Reynolds number 

"^Simulation time in orbits 

'^Surface density ratio of active to dead zones 

*^RMS y-component of force on protoplanet over orbits 25-100 

^Torque correlation time in orbits 



Table 2. Skew for 10 orbit intervals from 25 — 95torb for 64 x 256^ runs. Columns are 

labeled by the interval center time in orbits 



Rcm 30 40 50 60 70 80 90 

oo -0.283 0.463 0.076 -0.157 0.209 -0.026 0.028 

100 0.175 -0.158 -0.505 -0.113 -0.170 -0.651 -0.163 

30 -0.240 0.305 0.277 -0.123 0.254 -0.680 0.455 

3 0.097 0.278 0.673 0.145 0.048 -0.001 0.076 
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Table 3. Kurtosis for 10 orbit intervals from 25 — 95torb for 64 x 256^ runs. Columns are 

labeled by the interval center time in orbits. 





30 


40 


50 


60 


70 


80 


90 


oo 


0.657 


1.446 


0.460 


0.379 


0.555 


1.025 


1.274 


100 


0.910 


0.304 


0.425 


0.064 


0.675 


0.944 


0.439 


30 


2.066 


1.133 


1.057 


0.201 


0.573 


3.112 


1.223 


3 


0.362 


0.746 


2.883 


0.923 


0.921 


1.386 


0.258 



Table 4. Ratios of RMS density fluctuations and magnetic pressure. 



Run Name 




0.5 a 


{SPLgY-'/cl ^ 




32Rinf 


3.80( 


-2) 


1.16(-2) 


3.27 


32R100 


2.51( 


-2) 


5.93(-3) 


4.24 


32R30 


2.80( 


-2) 


2.80(-3) 


9.98 


32R3 


6.84( 


-3) 


3.12(-4) 


21.9 


64Rinf 


4.27( 


-2) 


1.93(-2) 


2.22 


64R100 


2.32( 


-2) 


4.01(-3) 


5.78 


64R30 


1.74( 


-2) 


2.28(-3) 


7.63 


64R3 


7.58( 


-3) 


2.22(-4) 


34.2 


128R30 


1.64( 


-2) 


1.26(-3) 


13.0 



'^RMS density fluctuations aX t — 50^0^6 

^variance of magnetic pressure normalized by aX, t — hQtorb 
'^ratio of cols a and b 



- 25 - 




M m m m m 



Fig. 1. — Space-time diagrams of magnetic energy averaged over horizontal (x — y) planes 
for runs with (from bottom to top) Rcm = oo, 100, 30, and 3. The dead zone is apparent 
in all models with finite Rcm- Note that the color scale is shifted to higher values for the 
bottom panel 
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Fig. 2. — Vertical profiles of Reynolds and Maxwell stresses for medium resolution runs 
with four different values of magnetic Reynolds numbers Rcm- Note that the Reynolds 
stress remains nearly constant across the dead zone, even as the dead zone depth increases. 




Fig. 3. — A resolution study for a Rcm = 30 model. Resolutions corresponding to dx = 
3.125 X 10^^ 1.563 x IQ-^, and 7.813 x 10"^ are labeled by = 32, 64, and 128. Variables 
shown include (a) the temporal RMS azimuthal gravitational force {Fy)f^ [stars) and the 
active-to-dead zone column density ratio J^a/^d (diamonds) versus resolution; and (b) volume 
averages of the effective dimensionless viscosity a = {{pUxUy) — {B^By)) / Pq (top), RMS Mach 
number (middle), and magnetic energy normalized to initial midplane pressure Pq (bottom), 
function of time in orbital units. 



-28- 




Fig. 4. — Time series of azimuthal gravitational force Fy in 64 x 256^ resolution runs. The 
force is scaled in units of 27rG'E, the vertical restoring force near the midplane. 
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Fig. 5. — RMS (left) and mean values (right) of the torque Fy for seven periods each of 
length lOtorb, labeled by the interval midpoint. The lOtorb means show weak fluctuations 
around the global time standard deviation, as much as a factor of ~ 0.2 in the Rcm = 3 
case, and a factor of a few in all cases. 
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Fig. 6. — Distributions of torque exerted on protoplanets in models with various i?eM- Each 
panel shows the distributions for seven periods, each of length lOtorb (from 25torb to OStorb), 
labeled by the interval midpoint. Note the decreasing range along the x-axis from model to 
model, indicating a narrower range of torque amplitudes for larger dead zone models. The 
distributions appear time stationary after the initial 10 orbits. 
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Fig. 7. — Power spectra of torque fluctuations Fy in four medium resolution models with 
increasing dead zone width. Note the decreasing value of the vertical scale from model to 
model. Very little power is seen at low frequencies. 
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Fig. 8. — Running time average of the torque Fy in the four medium resolution models with 
increasing dead zone size, starting at t = 25torb to eliminate effects of the initial transient. 
The dotted hne marks the zero point. The ordinate is a factor of ten smaller for the Rcm — 3 
case than for the others. In all cases, the running average approaches zero at late times. 
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Fig. 9. — Autocorrelation functions of the torque {Fy) fluctuations in the four medium 
resolution models with increasing dead zone width. All models have similar correlation 
times, roughly defined as the point where the curve first crosses zero with positive first 
derivative. 
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Fig. 10. — Cumulative torque contributions (5Fj, from spherical shells as a function of radial 
distance r from the box center shown at four random times in medium resolution models 
with increasing dead zone size (from top to bottom, Rcm = C)0, 100, 30, 3). In all models, 
the torque vanishes at or near r H, the disk scale height. The thin lines mark the zero 
point on the ordinate and r = H on the abscissa. 
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Fig. 11. — Turbulent overdensities 6p = p—po, where po is the initial hydrostatic distribution, 
in three orthogonal planes through the origin in each model with increasing dead zone width 
and normalized by the initial midplane density, Po(0). xz and yx cuts surround the yz 
cut. All sub-panels share a common color scale, given by the colorbar over the yz images. 
Clockwise from top left. Re = oo, 100, 30, 3. As the dead zone width increases, a clear 
transition from turbulent to wave-like behavior is visible at the midplane. 
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Fig. 12. — Derived diffusion parameter e estimated from turbulent torques versus dead zone 
to active zone column density ratio S^i/Eyi in all runs. Each resolution is labeled by its 
N,j. = 32, 64, 128, corresponding to our low, medium, and high resolution runs. For the low 
and medium resolutions, there are several Rcm values; for the = 128 case there is only 
one {Rcm = 30), marked by a diamond. The circled points correspond to the ReM — 30 
case at each resolution, allowing a resolution study for D versus E^/S^. The dead zone size 
should be more resolution dependent towards the right of the figure, as the active layers get 
thinner and thus increasing numbers of unstable modes are inadequately resolved. 



